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REDUCING  ALIASING  IN  THE  WIGNER  DISTRIBUTION 
USING  IMPLICIT  SPLINE  INTERPOLATION 
by 

G.F.  Boudreaux-Bar t el s  and  T.  W.  Parks 

I*  Introduction 

The  Wigner  Distribution  (WD), 

°° 

( t , to )  =  J  x ( t  +  t/2)x  (t  -  -r/2)e  (1) 

—CD 

of  a  continuous  time  signal  x(t)  is  a  bilinear  signal  transformation 
that  is  useful  for  analyzing  signals  whose  frequency  content  changes 
with  time.  The  WD  originated  in  the  field  of  quantum  mechanics,  and  has 
been  used  to  characterize  the  time— varying  signals  found  in  optical  com¬ 
munications,  speech  waveforms  and  radar/ sonar  patterns.  The  ambiguity 
function  of  x(t)  is  related  to  W^(t,w)  via  two  Fourier  transforms. 

Wigner  distributions  exhibit  several  desirable  properties:  1)  the  auto 
WD  of  a  real  or  complex  signal  x(t)  is  real,  2)  W^(t,co)  has  the  same 

time  and  frequency  support  as  x(t),  and  3)  shifts  in  time  or  frequency 
on  x(t)  produce  corresponding  shifts  in  W^(t,o>).  In  addition,  the  aver¬ 
age  frequency  of  W^(t,u>),  i.e.  fixed  t^,  is  equal  to  the  instantaneous 
frequency,  whereas  the  average  time  of  W^(t,(i)^)  is  equal  to  the  group 
delay  of  the  signal.  Integration  of  W^(t,(o)  with  respect  to  the  fre¬ 
quency  (time)  variable  for  a  fixed  time  t^  (frequency  co^)  yields  the 
instantaneous  power  (energy  density)  of  the  signal  at 


integral  over  the  entire  (t,w)  plane  is  equal  to  the  signal  energy. 
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W  (t,io)  can  be  viewed  as  the  Fourier  Transform  of 
x 

=  x(t  +  t/2)x  (t  “  t/2)  ,  (2) 

the  product  of  appropriately  shifted  versions  of  x(t)  and  its  complex 

* 

conjugate,  x  (t).  When  the  Fourier  integral  cannot  be  calculated  in 
closed  form,  digital  signal  processing  techniques  are  often  invoked  to 
obtain  numerical  approximations.  A  typical  signal  processing  approach 
is  to  sample  either  g^(t,x)  or  x(t),  window  the  resulting  sequence,  and 

then  compute  its  discrete  Fourier  Transform  (DFT).  However,  x(t)  must 
be  sampled  at  twice  the  Nyqui st  rate  [2]  to  avoid  aliasing  errors  when 
calculating  W^(t,w).  In  many  signal  processing  situations  where  under- 

sampling  has  occurred,  anti-aliasing  digital  filters  can  be  designed  to 
minimize  aliasing  errors  [4,5].  An  alternate  approach  is  to  view  the 
sequence  of  samples  as  part  of  a  spline  function  approximation  to  the 
continuous- time  signal  x(t).  The  Fourier  integral  of  a  spline  approxi¬ 
mation  can  be  easily  calculated  using  certain  predetermined  weighting 
functions,  called  diminishing  factors  [6],  to  reduce  aliasing  errors. 
In  many  cases,  the  spline  function  analysis  method  produces  a  smaller 
approximation  error  to  the  Fourier  integral  than  the  DFT  of  the  sequence 
of  samples  does. 

It  is  the  purpose  of  this  research  to  extend  the  theory  of  implicit 
spline  approximations  of  the  Fourier  integral  to  similar  approximations 
of  the  Wigner  distribution.  We  determine  the  two-dimensional  general¬ 
ized  diminishing  operators  needed  to  correct  the  aliasing  errors  gen¬ 
erated  when  only  an  undersampled  version  of  x(t)  is  available  to  approx¬ 
imate  W^t,^).  This  implicit  interpolation  process  "fills-in"  the 
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sampled  signal,  x(nT),  between  the  signal  samples,  increases  the  effec¬ 
tive  sampling  rate,  and  hence  reduces  the  effects  of  aliasing  when  cal¬ 
culating  the  Wigner  distribution.  The  spline  approximation  can  be  cal¬ 
culated  very  efficiently  if  one  takes  advantage  of  the  special  structure 
and  symmetry  of  Wigner  distributions. 

II.  Theoretical  Developments 

Classen  et  al  [2]  and  Chan  [3]  have  studied  the  effects  of  sampling 
and  windowing  x(t)  in  order  to  approximate  W^(t,o>).  If  x(t)  is  bandlim- 

i  ted  to  to  ,  i.  e. 
s 


X(o>) 


CO 

f  x(t)e^COtdt  =  0,  I  to  I  >o) 

J  —  s 

—CO 


2  jrf  , 
s 


(3) 


then  W  (t,to)  is  also  bandlimited  to  o>  .  The  DFT  can  be  used  to  calcu- 
x  s 

lated  (3)  on  a  discrete  grid  of  frequency  points  if  x(t)  is  finite 

length  and  sampled  every  T  ^  ~r~  seconds.  However,  twice  as  many  sam- 

s 

pies  of  x(t),  is  needed  to  avoid  aliasing  errors  when  calculating 
Wx(t,o>)  on  the  same  frequency  grid,  i.e.  T 

s 

If  the  signal  x(t)  has  been  sampled  at  the  Nyquist  rate,  then  one 
can  avoid  aliasing  errors  in  approximating  W^(t,w)  by  interpolating  the 

sampled  sequence.  However,  in  addition  to  the  computational  expense 
involved  in  interpolation  algorithms,  doubling  the  number  of  samples 
more  than  doubles  the  number  of  calculations  needed  to  approximate 
W^(t,w).  An  alternate  approach  to  minimizing  the  effects  of  aliasing  is 


to  view  the  undersampled  sequence  as  the  equally  spaced  knots  of  a 


5 


spline  approximation  to  the  continuous-time  signal  x(t). 

Let  the  discrete  signal  (x(n)}  represent  the  uniformly  spaced  sam¬ 
ples  of  the  continuous  time  signal,  x(t).  That  is, 

x(n)  =  x(nT)  ,  n  =  *..,  -2,  -1,  0,  1,  2,  ...  .  ^ 

The  Fourier  transform  of  this  sampled  sequence  is 

X(f)  =  Y  x(n)e"j2nfnT  .  (5) 

n  =  —co 

A 

Furthermore,  assume  x(t)  is  the  first-order  spline  approximation  of  x(t) 
that  interpolates  the  samples,  x(nT).  The  Fourier  integral  of  the 
spline  approximation 


X(f)  =  J  x(t)e  j2nftdt  =  X(f)  (6) 

—  CO 

is  easy  to  compute  [6]  since 

i(f)  =  X(f)  [-innf^T]2  (7) 

is  the  product  of  the  Fourier  transform  (5)  of  the  signal  samples,  and 

.  „  2 
the  "diminishing  factor"  [  ( sinnfT) /nfT]  .  This  diminishing  factor  is 

predetermined  according  to  the  order  of  the  spline  approximation  and 

weights  the  periodic  spectrum  of  the  undersampled  sequence,  x(n),  to 
produce  the  unaliased  spectrum  of  the  continuous  time  spline  approxima— 

A 

tion,  x(t).  This  technique  is  said  to  be  "implicit"  since  the  spline 

A 

function,  x(t)  is  never  actually  computed.  Note  that  since 
2 

[( sinnfT) /nfT]  is  the  Fourier  transform  of  the  triangular  pulse 
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p(t) 


Itl  <  T 
otherwi se 


(8) 


the  spline,  or  piecewise  linear,  approximation,  x(t),  can  be  written  as 


the  convolution  of  the  discrete  sequence  (x(n)}  and  the  triangular 
interpolating  function,  p(t): 


09 

x(t)  =  £  x(n)p(t-  nT) 

n=— <» 


J  x(nT)p(t  -  nT) 

n=>-co 


(9) 


In  order  to  apply  these  results  to  ffigner  distributions,  let 


WA(t,w) 


J  £(t  +  t/2)£  (t  -  t/2) e_Jtl)Tdr  =  W^t.co) 


(10) 


be  the  WD  of  the  spline  approximation,  x(t)*  Plugging  (9)  into  (10)  we 
obtain 


WA(t,w)  =  J  £  x(iT)x*(mT)W  (t  -  T,(o)  (11) 

x  i=— 00  00  ^ 

where  W  (t,w)  is  the  WD  of  p(t)* 

P 

Evaluating  (11)  at  t  =  nT,  we  find 


a>  a? 


WA(nT.u)  =  £  J  x[(n  +  i)T]x*[(n  -  m)T]e  T,«]  (12) 


m=-oo 


Since  W  tt,w]  =  0,  V  Itl  >  T,  (12)  reduces  to 
P  ~  “ 
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1 

WA(nT,w)  =  J 
x  q=-l 


{  ]jj  x(n  +  i)x((n 

i=-o» 


v  . »  —  j  <jo2  i  T- 
q)  -  l)e  ]x 


{e  j<J)qTW  [ qT/ 2 ,  u>] } 
P 


(13) 


Equation  (13)  is  similar  in  form  to  (7)  since  the  expression  in  the 
first  set  of  brackets  is  merely  the  discrete  cross  Wigner  distribution 

[2]  of  the  undersampled  sequences  x(n)  and  x(n  -  q) .  Equation  (13) 

could  be  calculated  by  appropriately  weighting,  or  diminishing,  three 

cross  Wigner  distributions*  However,  a  much  more  efficient  method  can 

be  obtained  by  taking  advantage  of  the  special  symmetry  that  results 

when  one  computes  the  DFT,  with  respect  to  n,  of  W^(nT,  g>). 

x 

III*  Implements  tion 

Define  the  convolution  sequence  (for  fixed  w^) 


h  (n) 

“o 


V  ~  VT  ~ 

J  (x(i)e  U  ) (x(n 


-jw  (n  -  i)T 
i)  e  °  ) 


CO 

=  ^  x(i)x(n 

i=>-oo 


-j(°02iT  j(o0nT 
i)  e  e 


Note  that. 


(14) 


h  (2n)  =  W  (nT, wA) 


(15) 


is  the  discrete  WD  of  the  samples  x(n)  evaluated  at  the  frequency. 


0 


Let 
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H  (z)  =  5  h  (n)z 

wA  UA 

0  n=-<*>  0 


(16) 


be  the  z-transfonn  of  h^  (n).  Combining  (13-14),  we  obtain 


WA(nT,<j>) 

x 


1 

h  (2n  -  q) W  [q  T/2,<i)]  . 
“  P 


If  we  define  the  z-transform 


(17) 


f(z,io0)  =  Y  WA(nT,w0)z  °  (18) 

n=-®  x 

and  insert  (16-17)  into  (18),  we  obtain 

5  {  5  H  (Z1/Vj’,i>2q/Vj”l5)»  [q  T/2.»J  (19) 

0  2  q~l  iio  “0  "  0 

Furthermore,  if  we  assume  that  x(n)  is  finite  length,  N,  and  also  assume 

(without  loss  of  generality)  that  x(n)  is  causal,  then  the  DFT  can  be 
used  to  evaluate  (16)  on  a  set  of  2N  discrete  normalized  frequency 
points,  i.e. 


j|gk  2S-1 

H  (e  2N  )  =  J 


n=0 


h  (n) 
“0 


•In 

J2N 


kn 


~  k  “o  -k  “o 

X(2NT  +  2n)X  (2NT  +  2nJ 


(20) 


where  X(f)  is  the  length  2N  DFT  of  the  nndersampled  sequence,  x(n)  and 


2  7T  2  7T 

Wq  is  one  of  the  discrete  frequency  points  =  0,  i  2NT*  i  2NT^ 1 


9 


2tt 

*  2Ni  '  •  Similarly,  (19)  can  be  computed  on  a  NxN  grid  of  normal¬ 

ized  frequency  points  by  noting  that 


“Mo 


{(£(- 


+  2m  — 


2NT 


2m  -  Ni 


2  NT 


)} 


* 

i  {Wp(0,j^m)  +  2(-l)icos(|gk)Wp[T/2,^|m]}  (21) 


since  W  (t,(o)  =  W  (-t,w)  for  p(t)  real  and 
P  P 

2N  DFT  of  x(n)  and  the  length  N  vectors 
only  be  computed  once  and  stored  in  memory 


symmetric.  Thus,  the  length 

( 0 , ~m )  and  W^[T/2,j~“m]  need 
.  Then  for  each  m  =  0,  1, 


.2 n 

th  J2nk  2n 

....  N  -  1,  the  m  column  t(e  ,^pn) ,  k  =  0,  1,  ....  N  -  1  and  the 

length  N  inverse  DFT, 


.  2n  .  1  Nf1  . ,  ^2Nk  2n  j  N  kn 

ff (nT' NT™  =  N  JQ  #(e  '^)e 


(22) 


are  computed  using  Fast  Fourier  Transform  (FFT)  techniques  to  effi¬ 
ciently  obtain  the  spline  approximation  of  the  Wigner  Distribution. 


TV.  Examples 


The  length  N  =  41  equal-ripple  finite  impulse  response  in  Figure  1 
was  designed  using  the  Remez  Exchange  algorithm  [7].  Its  frequency 
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response,  shown  in  Figure  2  was  designed  with  a  unity  passband  from 

Fp  *=  0.0  to  Fp  =  0.15.  The  stopband  extends  from  normalized  frequency 
1  2 

Fe  =  0.2  to  Fc  =  0.5.  The  passband  and  stopband  errors  were  equally 
bl  b2 

weighted.  Note  that  the  impulse  response  in  Figure  1  was  sampled  more 
than  twice  as  fast  as  required  by  the  Nyquist  criterion  in  order  to 
avoid  aliasing  errors  when  calculating  its  discrete  ffigner  distribution 
shown  in  Figure  3.  The  WD  shown  in  Figure  4  is  that  of  the  impulse 
response  in  Figure  1  decimated  by  a  factor  of  2.  Even  though  the 
reduced  sampling  rate  satisfies  the  Nyquist  criterion,  it  is  not  fast 

enough,  i*e.  it  does  not  meet  the  requirement  that  T  ^  ,  to  avoid 

s 

aliasing  errors  in  the  WD.  Finally,  the  WD  of  the  spline  approximation 
of  the  decimated  impulse  response  is  shown  in  Figure  5.  Note  that  the 
effects  of  aliasing  have  been  removed.  The  first  order  spline,  or 
piecew i se-linear,  approximation  did  very  well  approximating  the  low  fre¬ 
quencies  and  only  moderately  well  approximating  the  high  frequencies,  as 
would  be  expected.  Higher  order  splines  can  be  used  to  obtain  better 
approxima  tions. 

V.  Conclusions 

Implicit  spline  interpolation  of  a  sampled  sequence  has  been  used 
to  approximate  the  Wigner  distribution  and  reduce  the  effects  of  alias¬ 
ing  errors  that  occur  when  a  signal  is  under sampl ed.  Recall  that  sig¬ 
nals  must  normally  be  sampled  at  twice  the  Nyquist  rate  in  order  to 
avoid  aliasing  when  calculating  its  WD.  The  implicit  approximation 
technique  is  an  efficient  alternative  to  oversampling  or  explicitly 
interpolating  a  given  sampled  signal.  The  symmetries  resulting  in  the 
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double  frequency  plane  in  (21)  can  be  utilized  to  compute  both  the 
discrete  Wigner  distribution  and  the  spline  approximation  of  the  WD  very 


efficiently 
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Figure  Captions 


Figure  1. 

Impulse  response  of  a  length  N  =  41  equal-ripple  lowpass 

filter  designed  with  unity  passband  from 

F0  =  0.20  to  Fd  =  0.15  and  stopband  Fc  =  0.20  and 

F1  P2  bl 

Fc  =  0.5.  The  sampling  interval,  T,  is 
b2 

assumed  to  be  equal  to  1.0. 

Figure  2. 

Frequency  response  of  a  length  N  =  41  equal-ripple 

lowpass  filter  designed  with  unity  passband  from 

Fp  =  0.0  to  Fp  =  0.15  and 

1  2 

stopband  ranging  from  Fc  =  0.20  to  Fc  -  0.5. 

bl  S2 

Figure  3. 

Time  (n)  vs.  frequency  (u>)  plot  of  the  Wigner 

distribution  (WD)  of  the  impulse  response  shown  in 

Figure  1.  (Sampling  interval  T  =  1.0) 

Figure  4. 

Time  (n)  vs.  frequency  (to)  plot  of  the  WD  of  the 

decimated  lowpass  impulse  response.  (Sampling  interval 

T  =  2.0) 

Figure  5. 

Time  (n)  vs.  frequency  (w)  plot  of  the  WD  of  the 

first-order  spline  approximation  using  the  decimated 
lowpass  impulse  response  sequence  as  the  equally  spaced 
spline  knots.  (Sampling  interval  T  =  2.0). 
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FIGURE  3 
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FIGURE  5 


